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Technical Field 

This invention concerns a method for estimating the frequency of a 
single frequency complex exponential tone in additive Gaussian noise. In 
another aspect the invention is a frequency estimation program for estimating 
the frequency of a single frequency complex exponential tone In additive 
Gaussian noise. In further aspects, the invention is computer hardware 
programmed to perform the method. 

Background Art 

Earlier work to estimate the frequency of a single frequency complex 
exponential tone in additive Gaussian noise uses the fast Fourier transform 
(FFT) algorithm. The initial work on this topic was introduced by Rife and 
Boorstyn [1-3]. This paper introduces an algorithm employing the FFT, which 
produces an estimate of the frequency with extremely low variance of the error. 
The variance of the frequency estimate is independent of the frequency of the 
signal. The algorithm has a low computational complexity implementation. 

The received signal, r[n], is given by, 

r[n] - s[n] + n[n], for n = 0,1 ,2 N-1 (1 ) 

where: 

s[n] = Ae^ T «, 

{V Mlo ' 1 is a set of independent, complex, zero mean, Gaussian random 

variables with variance a 2 , 
7 R [n] = Re{7[n]}, 

7 I [n] = Imag{7[n]}, 

f is the frequency of the tone, 

T s is the sampling period, 



<x 2 

— = var[77 R [n]] =var[7,[n]] 
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and, A is the signal amplitude. 

The sampling frequency, f 8l is given by, 

f s =~ samples/s (2) 

*S 

The signal to noise ratio of each complex signal plus noise sample is given by, 
A 2 



SNR=- 



(3) 



Rife and Boorstyn [1-4] suggest a method of estimating f by using a FFT. 
It is assumed that 0 < f < f s . First, a coarse search is performed. Under 
noiseless conditions, the absolute value of the FFT output coefficient 
corresponding to the bin centre frequency closest to f will be maximum over the 
set of absolute values of the FFT output coefficients. The coarse search, 
performed by the FFT, narrows the frequency uncertainty, to 

^ Hz, where an N point FFT is used. Then, a fine search method is used to 

further reduce the frequency uncertainty. A secant method is used to compute 
the estimate of f by successful approximates. 

Define, 



(4) 





"r(0) 




"Y(0) 




KD 




Y(1) 


r = 


r(2) 


, Y = 


Y(2) 




_r(N-1) 




Y(N-1) 



where, Y =jRFT(r) and FFT(.) is the Fast Fourier Transform Operator. 

Then the Rife and Boorstyn coarse search is, 
k nmx =max- 1 {|Y(k)|:0£k£N-l} . 



(5) 



and, 
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(6) 



where, 

f 0 is the coarse frequency estimate in Hz. 

Numerous other frequency estimation approaches have been suggested 
in the literature [5 -1 0]. 

Disclosure of the Invention 

In a first aspect the invention is a method for estimating the frequency of 
a single frequency complex exponential tone in additive Gaussian noise, 
comprising the steps of: 

performing the fast Fourier transform (FFT) on the tone; 

estimating the frequency as the frequency corresponding to the largest 
FFT output coefficient magnitude; 

computing a discriminant which is proportional to the frequency error in 
the initial frequency estimate using modified coefficients of the discrete Fourier 
transform (DFT) with center frequencies plus one half and minus one half of the 
FFT bin spacing relative to the initial frequency estimate; 

mapping the value of the discriminant into the estimate of the frequency 
error in the initial frequency estimate using a mathematically derived function; 

adding the estimate of the frequency error to the initial frequency 
estimate to get a first interpolated frequency estimate; 

computing a further discriminant which is proportional to the frequency 
error in the first interpolated frequency estimate using modified coefficients of 
the discrete Fourier transform (DFT) with center frequencies plus one half and 
minus one half of the FFT bin spacing relative to the first interpolated frequency 
estimate; 

mapping the value of the further discriminant into the estimate of the 
frequency error in the first interpolated frequency estimate using the 
mathematically derived function; and 

adding the estimate of the frequency error in the first interpolated 
frequency estimate to the first interpolated frequency estimate to get a second 
interpolated frequency estimate. 
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The first interpolated frequency estimate is quite accurate because it is in 
a region of relatively low noise induced frequency error. The method generates 
an unbiased, low error variance estimate of the frequency. The performance of 
the method, above the signal to noise ratio threshold, is about 0.06 dB above 
5 the Cramer-Rao lower bound. The method is ideally suited to be utilised in a 
number of communications, signal processing and biomedical applications. 
The method is easily implemented in hardware or software with low 
computational overhead. 

In theory, this technique of iteratively deriving an interpolated frequency 
10 estimate and then, using the frequency discriminant, a more precise frequency 
estimate can be continued infinitely times until a fixed point (or solution) occurs. 
At this fixed point, the discriminant function has zero value. 

Several functions have been identified to compute the discriminant. In 
practice, different functions may require a different number of iterations to 
15 essentially converge to a fixed-point solution. However, discriminant functions 
defined by a wide class of functions using two DFT coefficients as the input 
converge to the same solution and therefore exhibit identical noise 
performance. 

A first example of the discriminant, or distance metric, of frequency 
20 estimation error is: 

where, e = fT s — =s- (1 0) 

N 

and, 

25 s = fr s -^!!« 

N 

So for the initial frequency estimate using the FFT, f 0 T s = kjS53L and e = 0. 

N 

Other examples of the discriminant having the properties required for the 
algorithm include: 



30 D = ii^-i«r >forr>0>f 

r\J3\'+\ar' 7 



Y\P\ r +\a\> 
and in particular, 
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D _i |/?! 2 -l<*l 2 

2|/?| 2 +M 2 



and 



D = Re[- 



J3 + a 



r] 



where Re[J is the real part and * denotes the complex conjugate. 

In these equations, p and a are the modified DFT coefficients defined by, 



It is also possible to define discriminant functions which use more than 
two DFT coefficients to obtain further improvements in frequency estimation 
performance in additive Gaussian noise relative to discriminants that use only 
two DFT coefficients. An example where 2M+2 coefficients are used, where 

0<M<,—-\ and the FFT coefficients are used in the discriminant with 

optimal weighting coefficients obtained by using the concept of matched 
filtering is, 




and, a = Y(k max -I) = 2r(n)e 



D=Re[ 




m)modN]} 



] 



£ {C nc 1,^^ +- + m)modN] + Y*[(k max -±- 

m=0 L o" mJ ^ Z 



m)modN]} 



where, 0£M z 1, mod N indicates modulo N, 

2 

and, where, * denotes complex conjugate. 
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and, Y(k max +^+m) and YOc^ -^-m) are the modified DFT coefficients given 
by. 



i±2 -I2fm i 



Y(k max+ i + m) = £r(n)e" J2 ' m ~^ 

Z n=0 

5 and, YCk^ - i - m) = ^m^ m ~^~ 

The discriminant using more than two DFT coefficients may be used in 
the last iteration to obtain additional frequency accuracy. In a similar manner, 
discriminant functions may be formulated which use more than two DFT 
10 coefficients and less or equal to all N FFT coefficients. 

Additional frequency accuracy may be obtained by computing the 
frequency discriminant recursively until convergence for the frequency estimate 
is reached. 

Convergence for the frequency estimate may be reached after zero to 
15 three iterations, depending upon the specific discriminant used and the signal 
to noise ratio. 

In any iteration, the frequency discriminant may be computed using any 
one of the functional forms: 

y may vary on each iteration. 

1 0 — ce * 

25 A f m (r) = — Re[ m » ] f s where| Re[ ] deno tes the real part and ' 

denotes the complex conjugate 

In general, the frequency incremental shift Af m (r) is related to the 
previously defined frequency discriminant, D, by, 

Af m (r)^D 
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The frequency discriminant may be driven to zero input and output 
values by either modifying the frequency of the DFT coefficients or frequency 
translating the signal. Signal frequency translation may be achieved by 
multiplication of the signal by a locally generated complex exponential signal. 
5 The advantage of frequency multiplication of the signal is that the algorithm 
may be implemented with a standard hardware, software, or combination 
hardware/software FFT. This FFT may be highly optimized for one or a 
multiplicity of processors operating as a system. 

The process for obtaining additional frequency accuracy may be scaled 
10 to save multiplies by scaling the frequency estimate during recursion. The 
process may involve a final step of multiplying the scaled frequency estimate 
fm+iT, with the sampling frequency f s to remove the scaling from the frequency 

estimate. 

In a second aspect, the invention is a frequency estimation program for 
15 estimating the frequency of a single frequency complex exponential tone in 
additive Gaussian noise, wherein the frequency estimation program has 
functionality to perform the method. 

In a third aspect, the invention is computer hardware programmed to 
perform the method. The hardware may comprise a DSP processor chip, or 
20 any other programmed hardware. 

Brief Description of the Drawings 

Examples of the invention will now be described with reference to the 
accompanying drawings, in which: 
25 Fig. 1 is a graph that illustrates the FFT Coefficients, where the signal 

frequency is closer to the lower FFT frequency than the higher FFT frequency; 

Fig. 2 is a graph that illustrates the FFT Coefficients, there are two equal 
peak coefficients and the signal frequency is halfway between; 

Fig. 3 is a graph that illustrates the FFT Coefficients, where the signal 
30 frequency is closer to the upper FFT frequency than the lower FFT frequency; 
Fig. 4 is a Flow Diagram for the Frequency Determination Algorithm; 
Fig. 5 is a graph that illustrates the ratio of the variance of the normalized 
frequency error, s-i, to Cramer-Rao Bound variance in dB as a function of 
the FFT length, N; and 
35 Fig. 6 is a graph that illustrates the variance of the normalised estimator 

frequency error estimate against the frequency error for the first interpolation. 
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Simulations of the invention show the rms frequency error performance 
of the algorithm vs SNR in dB, for different values of N. 

Figures 7-12 include curves for one interpolation, two interpolations, and 
the Cramer-Rao Bound, where: 

Fig. 7 is a graph showing RMS normalised frequency error vs SNR for 

N=2; 

Fig. 8 is a graph showing RMS normalised frequency error vs SNR for 

N=4; 

Fig. 9 is a graph showing RMS normalised frequency error vs SNR for 

N=16; 

Fig. 10 is a graph showing RMS normalised frequency error vs SNR for 

N=64; 

Fig. 1 1 is a graph showing RMS normalised frequency error vs SNR for 
N=256; and 

Fig. 12 is a graph showing RMS normalised frequency error vs SNR for 
N=1024. 

Fig. 13 is a Flow Diagram for the Frequency Determination Algorithm 
using a fixed number of iterations stopping rule. 

Fig. 14 is a Flow Diagram for the Frequency Determination Algorithm 
using a magnitude of the frequency error discriminant stopping rule. 

The two-interpolation case essentially achieves the performance of the 
infinite interpolation case. 

Best Modes of the Invention 

Referring first to Figs. 1 to 4, the received signal r[n] is given by: 
r[n] = s[n] + T)[n], for n =0,1,2, .... N-1 ' (-j) 

s[n] = Ae j2,rfi, V 
where: 

torn]}? 4 

is a set of independent, complex, zero mean, Gaussian random variables 

%[n] = Imag{i7[n]}, 

f is the frequency of the tone, 

T s is the sampling period, 
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and, A is the signal amplitude. 

A fast Fourier transform is performed and the sampling frequency, f Sl is given 
5 by, 



Then, an initial frequency estimate f 0 is taken as the frequency 

corresponding to the largest FFT output coefficient magnitude. A discriminant 
which is proportional to the frequency error in the initial frequency estimate f 0 
10 is computed using modified coefficients a 0 , fi 0 of the discrete Fourier 

transform (DFT) with center frequencies plus one half and minus one half of the 
FFT bin spacing relative to the initial frequency estimate ^ . 

The value of the discriminant is then mapped into the estimate of the 
frequency error in the initial frequency estimate f 0 using a mathematically 
15 derived function. 

The estimate of the frequency error is added to the initial frequency 
estimate f 0 to get a next interpolated frequency estimate f t . 

The process is then repeated, using the next interpolated frequency 
estimate f, and computing a new frequency discriminant to produce a next, 
20 more precise, frequency estimate f 2 . 

THE FREQUENCY INTERPOLATION DISCRIMINANT 
It is assumed that the signal to noise ratio (SNR) is sufficiently high such 
that the largest magnitude FFT coefficient corresponds to a frequency closest 
25 to the signal frequency. This assumes that the signal to noise ratio is 
sufficiently high that the probability of the statistical outlier event of a noise only 
FFT bin magnitude being larger than a FFT bin containing both signal and 
noise is negligible. 




samples/s 



(2) 



30 



Define, 

N-l 



2N 



1 



)n] 



(7) 
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N-l 



P = 2 r[n]exp[-j2*(f + -Un] 



(8) 



Then the discriminant, or distance metric, of frequency estimation error is 
defined as, 

+ (9) 
where, s = fT 3 (10) 

and, 

k 

s N 

For the initial frequency estimate using the FFT, f 0 T s = ^m. and £ = 0. 

N 

in the noiseless case, 

-1, S—S — 

2N 

1 - 1 

1, e-s = — 

2N 

D(e,e) is a monotonically increasing function qU-e. Therefore, 
each D(e,i) , there is a unique inverse mapping to s-e. Clearly, D(e,i)) may 
be used as a discriminant for fine frequency interpolation between FFT bin 
centre frequencies. 

There exists some functional relationship such that, 

?,T.=^ +I *D(M)] . (12) 

where, \y( .) is a monotone increasing function. 

y/{.) is called the frequency interpolation function and 

f, is the first interpolated frequency estimate. 

The requirement that f, has zero error in the noiseless case is, 

y^D{e y £)] = e - e , for - 1 £ D £ 1 . Therefore, ys~ l (s-i) = D(e, i) . 

THE FREQUENCY INTERPOLATION FUNCTION 

k 1 k 1 

Assume that, — — ^ fr £ +— . 

N 2N N 2N 
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Tl»n.- 5 L*.-**-L. (13 ) 

r[n] = s[n] = Ae J2s,hT ' > n = 0,1,2,..., N-l (14) 

Without loss of generality, assume that £ = 0. Also assume the noise 
5 free case. 

The FFT output coefficients are given by, 

YCO'A ^ , k = 0,l,...,N-l 

_ Ae j«tN-o sinQreN) ( 1 5 ) 
sin(TZE) 



The discriminant can be expressed as , 
10 ^- |Y(k - + i )HY(k -4' lsin(^^)l-|sin(^-^), 
I Y(k max • + 1) | + 1 Y(k _ - 1) | | sin(^ + ^) | + 1 sin(^ - JL) | 

After some trigonometric simplification, 

tan(— ) 

D(s,0) = 2L. 

tan(a») 



(16) 



(17) 



15 This inverse mapping from T>(e,e) to e- e can be obtained as, 

7t 2N 



5-5=^[D(5^)]=ltan 1 [D(^,^)tan(^)] (18) 



Then the first interpolated frequency estimate, f„may be obtained, 
where, 

20 fiT3=^ +ff =^ + ^) = ^+^tan 1 [D( ff ,i)taa(^)] (19) 
V° "I V**^ (20) 



25 



The implication is that a two point (N=2) FFT is sufficient to obtain zero 
error frequency determination in the noiseless case. However, the Cramer-Rao 
lower bound is relatively large for N=2 and the SNR threshold is relatively large. 
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There is a motivation to use larger N to reduce the rms frequency estimation 
error. However, the implication of larger N is increased computational 
complexity and longer delay time to obtain the transform results. It is desirable 
to obtain low frequency estimation error with the smallest possible N. 

Define, 
^,(D) = ^- ,for-l£D£l. 



For large N or for any N and small |D|, the function yr x (D) closely 
10 approximates ^(D) . 

$T.- lim ^-+^(D) = lim *^+[ 2 

lYOc^^KlYOc^-I)! 2N 

N — >■ oo N -» eo (21) 



FREQUENCY ESTIMATION BY ITERATION 
For the case of r[n] consisting of signal plus noise, the noise will cause a 
15 perturbation of D(e,e), and some error in the frequency estimate will result 

Although an algorithm for the exact frequency determination in the noiseless 
case has been presented, it will be shown that the noise performance of this 
algorithm improves substantially when \s-e\ is close to zero. Since the 
discriminant T>(e,e) can be used to get an interpolated frequency estimate with 

20 less than one half FFT bin size error, it then follows that the algorithm can be 
iterated to move \ e-s\ towards zero and |D(*,£)| towards zero. In this way 

the variance of the frequency estimator output can be reduced. 
The iterative algorithm is defined below. 
Define, 

2r[n]exp[-j2^ + ^)] - 2r[n]exp[-j2^ - -^-)] 



25 N-i f -I N-l f 1 

£r[n]exp[-j2*(-f- + — )] + £r[n]exp[-j2^ - 



Define a monotone increasing function of 
D, <*(D), such that #0) ^ 0, <*(1 ) = ~,and #-1 ) = ~. 
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and, 



c _ -pT "•max 

N 



f =f T -^JS22L 

m m 8 N 

5 The iterative algorithm is defined by, 
k 

N 

which implies that £ 0 = 0 ■ 

10 Then, 

f/T s =f 0 T s +#D(*,s 0 ] 

■ • • 

f m T 8 =f m . 1 T s 4^[D(^,^ m . 1 ] 
15 and, 

m— >oo 
=lim(D m ) 
m->oo 

The steady state frequency estimate at the end of the iteration is f^. 
20 Define, s k =f k T a , for k=0,1 ,2,3 

Then, 

k 

A -p Tp max 

*- -I » *• ^~ 

and, the normalized frequency error is, 

4,-*=(L-f)T 8 

25 The iteration may be viewed as a convergence to a fixed point of the equation, 

K = g(.e n -\) = -<Z>[D(*>£m-i)]» for m £1 
1 ... 1 
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Theorem 1 below is referenced, [21]; 

Let g(x) be a continuous function on [ — — ,— ] and g([a,b])c [ — — — i 

2N 2N J/ 1 2N 2N 

5 Furthermore, assume that there is a constant 0 < X < 1, with, 

l#(x)-g(y)M|x-y|, forallx,ye[a,b], 

Then, 

x = g(x) has a unique solution x^ in [a,b], 
also, 

the iteration x„ =g(x n _, ),for n £ 1 will converge to 

x^for any choice of x„ e[a,2>], 

and, 

l x „-x n |<I— |x,-x 0 | 

In the situation under analysis, for fixed {T[n]}™ 1 and s, $>(.)is a function of i. 

10 

g(e) = 
and, 

|g(x)-g(y)N|x-y||l-^lM|. 

x-y 

If0< ^(x)-^(y) <2 
x-y 

Then, 

|g(x)-g(y)|<A x-y|, 

where,A = max|l- ^ x >-^> | for x,ye [-—,—] 

x-y ^ J 1 2N'2N J 

Using Theorem 1, it follows that the iteration will always converge to a fixed 
point under the appropriate conditions. 

15 Also, ^(DJ = 0,and from the properties of >(.), D ra = 0. 
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The fixed point solution, ^satisfies, 

1 2r[n]exp[-j2^ + -L)n] | - 1 2r[n]exp[-j2;r(^ - -L)n] | 

OO vt « A — A — 



| £r[n]exp[-j2;r(-f- + ~-)n] | + 1 f>[n]exp[-j2* - -L)n] | 
n=o t» 2N £3 f 8 2N 



The two previously defined functions y/(J>) and ^(D) fulfil the 
5 requirements of ^(D) and may be used in the iteration to obtain f a . While 
y(D) iteration will tend to converge more rapidly than ^,(D) iteration, both will 
yield identical values of f,. However, evaluation of ^(D)has lower 
computational complexity than evaluation of ty(D). There is performance 
advantage in using ^(D)when the computation is limited to a few iterations. 

10 

ADDITIONAL FREQUENCY ACCURACY 

Assuming that the SNR is sufficiently high, it is highly probable that 

f ~m tio + 2N ] ' 7,118 is an above thresh0,d condition [1]. A fine interpolation 
15 is obtained to improve the frequency accuracy. 

The normalized frequency estimate fr, is computed recursively in order 
to save computational complexity. 

After computation of fr 8 , f is obtained by (fr,)f 8 =f 

20 First Algorithm 

Referring to Figures 13 and 14, a first algorithm is provided to improve the 
accuracy of the frequency estimation. At step 1, the N point complex FFT is 
computed. The FFT output coeffficents are Y(k), 0<k^N-1 . 

25 At step 2, the peak search to find kmax is: k mat =max , [|Y[k]|:o^k^N-i]. 
At step 3, the initial frequency estimate is computed by: f 0 = k^f. 
At step 4, recursion is started at m=0. 

At step 5, the DFT coefficients for the m:th frequency estimate are computed: 

30 « m =|>]e ^ »> 

n=0 
n=0 
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The frequency discriminant is then computed: 
5 The m+1 :th frequency estimate, f m+ i , is computed: 

fm + i=f m + Af m (r) 

At step 6, convergence for the frequency estimate is reached if |f m+1 -f m | is 
sufficiently small. If there is convergence, then the frequency estimate is f m+1 . 
10 If convergence for the frequency estimate has not been reached, then m is 
incremented by 1 and step 5 is repeated. In practice, the algorithm will 
converge for m=2. Therefore, only f m for m=0, 1 , and 2 need to be computed 
(two iterations) which means that the frequency estimate is i 2 . 

15 Second Algorithm 

To handle large N, a modification to the first algorithm is made. The 
modification is made to the step of computing the frequency discriminant: 



20 




S 1 taa-lrHml-lg.nl (jL-\f 
* 1/U + IOnJ UNj^ 



s 1 r l>gm|-|«ml 1f 

2N L lAJ + |a m | J 8 



Third Algorithm 

A third algorithm is provided to improve the frequency accuracy. At step 1, the 
the N point complex FFT is computed. The FFT output coeffficents are Y(k), 

25 

At step 2, the peak search to find kmax is: k,^ =max- 1 [|YTk]|:o^k^N-i] 



1 

At step 4, the DFT coefficients for the initial frequency estimate are computed: 



At step 3, the initial frequency estimate is computed by: f 0 =^s.f 8 

N 



m -j*r»A_L) 



30 « 0 =2>Me Vf . 2N ' 



n=0 



n=0 
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The frequency discriminant is then computed: 

n IA>l + l«ol 2N 

5 The first interpolated frequency estimate Is computed: 

f,=f 0 +Af 0 (r) 



10 



At step 5, recursion is started at m=1 

At step 6, the DFT coefficients for the m:th frequency estimate are computed: 



:£r[n]e 2N 

D=0 



A»=2>[n]e ' 

n=0 



15 The frequency discriminant is then computed: 

Af m (r)= — rlflgf-lgml 2 -, t 
m(> 4NV m l 2 + |a m | 2] ^ 

The m+1 :th frequency estimate, f m+ i , is computed: 

f m+I =f m + Af m (r) 

20 

At step 7, convergence for the frequency estimate is reached if |f m+1 -f m | is 
sufficiently small. If convergence has been reached, then the frequency 
estimate is If convergence for the frequency estimate has not been 
reached, then m is incremented by 1 and step 6 is repeated. In practice, the 
25 algorithm will converge after f 2 is evaluated. Therefore, only f m for m=0, 1 , and 
2 need to be computed. This algorithm is less computationally complex than 
the first algorithm and has essentially the same convergence properties in the 
recursion. 

30 Fourth Algorithm 

A fourth algorithm is provided to improve the frequency accuracy. At step 1 , the 
N point complex FFT is computed. The FFT output coeffficents are Y(k), 

At step 2, the peak search to find Iw is: =max- 1 [|Yrk]|:0^k^N-i] 
35 At step 3, the initial frequency estimate is computed by: f 0 = kmax f s 

N 

At step 4, the DFT coefficients for the initial frequency estimate are computed: 
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h 1 



a 0 = £r[n]e '* 2N 

n=0 



N4 .J2irn(^.| U 
n=0 

The frequency discriminant is then computed: 

Af _ * r l>go l-I^Q | 
0 2nVoI + I«oI J " 

The first interpolated frequency estimate is computed: 

fi=f 0 +Af 0 

At step 5, recursion is started at m=1 

At step 6, the DFT coefficients for the m:th frequency estimate are computed: 

N-l -n„ n /im — L) 
« m =£r[n]e * 2N 

n=0 

N-l -j27rnf ^ 1 1 ) 

15 f ' 2N 

n=0 



10 



The frequency discriminant is then computed: 

Af _ 1 r IAnt 2 -|« m | 2 1f 

m 4NV m i 2 +i« m r s 



20 The m+1 :th frequency estimate, f m+1 , is computed: 

4«=f m + Af m (r) 

At step 7, convergence for the frequency estimate is reached if |f m+ i-f n | is 
sufficiently small. If convergence has been reached, then the frequency 

25 estimate is f m+1 . If convergence for the frequency estimate has not been 
reached, then m is incremented by 1 and step 6 Is repeated. In practice, the 
algorithm will converge after f 2 is evaluated. Therefore, only f m for m=0, 1, and 
2 need to be computed. This algorithm is less computationally complex than 
the first, second or third algorithms and has essentially the same convergence 

30 properties in the recursion. There is reduced computational complexity in the 
computation of Af,(r) because of the elimination of the need to compute square 
roots in the evaluation of the absolute value of complex variables. 
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Fifth Algorithm 

A fifth algorithm is provided to improve the frequency accuracy. At step 1 , the N 
point complex FFT is computed. The FFT output coeffficents are Y(k), Osk^N- 
1. 

5 At step 2, the peak search to find Wis: k max ^max 'tlYMIiO^k^N-i] 
At step 3, the initial frequency estimate is computed by: f 0 = kmax f 8 
At step 4, recursion is started at m=0 

At step 5, the DFT coefficients for the m:th frequency estimate are computed: 

The frequency discriminant is then computed: 
15 Af m (r) =^-[|^p~|^p3 f, , where . Ym is a set of non-negative constants. 

V m isaconstant,y m >0 

The m+1 :th frequency estimate, 4+, , is computed: 

W=4 + Af m (r) 

20 

At step 6, convergence for the frequency estimate is reached if |f m+1 -f m | is 
sufficiently small. If convergence has been reached, then the frequency 
estimate is f m+ , . If convergence for the frequency estimate has not been 
reached, then m is incremented by 1 and step 5 is repeated. 

25 

Scaled Computational Version of an algorithm 

The frequency scaled frequency estimate can be computed and then multiplied 
by f s . This process is described using the example of the first algorithm but can 
similarly be done for all the algorithms. The scaled computational version of the 
30 first algorithm is more computationally efficient as it saves multiplies. 

At step 1, the N point complex FFT is computed. The FFT output coeffficents 
areY(k),0<k<N-1. 

At step 2, the peak search to find Iwls: k,^ ^max-'EiYrkjlrO^k^N-i] 
35 At step 3, the initial frequency estimate is computed by: f 0 T s =»jag 

N 

At step 4, recursion is started at m=0 

At step 5, the DFT coefficients for the m:th frequency estimate are computed: 
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a m =2;r[n]e 2N 

n=0 

5 The frequency discriminant is then computed: 
The m+1:th frequency estimate, f m+1 , is computed: 

10 f m+ ,T,=f m T 8 +Af m (r)T 8 

At step 6, convergence for the frequency estimate is reached if |f m+1 -f m | is 
sufficiently small. If convergence has been reached then the scaled frequency 
estimate is f m+l T s . If convergence for the frequency estimate has not been 
15 reached, then m is incremented by 1 and step 5 is repeated. 

As a final operation, one additional multiply is require to remove the scaling 
from the frequency estimate. 

20 tfm+iT 8 )f, = f m+1 since T s f s =1 , where f m+ i Is the m+1 :th frequency estimate. 

Sixth Algorithm 

A sixth algorithm is provided to improve the frequency accuracy. The sixth 
algorithm uses any of the previously defined functional forms for A/ m (r) for any 

25 step. The difference between the sixth algorithm and the other algorithms types 
is that the frequencies of the two modified DFT coefficients are not changed. 
Instead, the centre frequency of the signal is modified by multiplying the 
previously defined signal by e ~ J 2nn *fm T s t0 obtain 
)x[n]e~ J2 * n * fmT ' }. The effect of this multiplication is to frequency 

30 translate the signal by - A/ m Hertz. The frequency discriminant is driven to zero 
input and output values by either modifying the frequency of the DFT 
coefficients or frequency translating the signal. The principle of driving the 
discriminant to zero by recursion is the same. The frequency error performance 
Of the algorithm as a function of signal to noise ratio is the same whether the 

35 signal is frequency translated or whether the DFT coefficients are frequency 
shifted. 
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There are situations where standard FFT or DFT functions are available in 
hardware, software, or combined hardware and software configurations. These 
FFT and DFT functions are highly optimised for their respective signal 
processors and are run at very high computational efficiency. Often, parallel 
5 processing for multiple processors is utilised extremely effectively. In these 
cases, the technique of frequency translation of the signal is of considerable 
implementation benefit. Very efficient computation is achievable. Frequently, an 
optimised large N point FFT runs faster on a parallel processor than the 
computation of two DFT coefficient 

10 

For the sixth algorithm, z 0 (n)= r(n\ n = 0,1,2,3,. . . is initialised. 

At step 1, the N point complex FFT is computed. The FFT output coefficients 
areY(k), 0SksN-1. 

At step 2, the peak search to find k^is: k,^ =max 1 [|Y[k]|:0^k^N-i] 
1 5 At step 3, the initial frequency estimate is computed by: f 0 - k ™» f s 
At step 4, recursion is started at m=0 N 
At step 5, the DFT coefficients for the m:th frequency estimate are computed- 

The frequency discriminant, A^(r), is then computed for any of the functional 

forms as a function of a m and fi m . 

The m+1:th frequency estimate, f m+1 , is computed: 

25 f m+1 =f m + Af m (r) 

At step 6, convergence for the frequency estimate is reached if |f m+1 -f m | is 
sufficiently small. If convergence has been reached, then the frequency 
estimate is f m+1 . If convergence for the frequency estimate has not been 
30 reached, then the signal is frequency translated (by complex time domain 
multiplication): 

^ + i(»)=^(«>" y2 ^ r ',/or « = 0,1,2,..., JV-1 
then m is incremented by 1 and step 5 is repeated. 

35 For convergence, all algorithms have the same steady state performance. 



WO 2004/005945 PCT/AU2003/000862 

22 

In the noiseless case, the first algorithm gives the exact frequency for one 
iteration, which Is of great benefit in some applications. The exact frequency is 
obtained for any N, including N=2. The reason is that an exact functional 
mapping from the magnitudes of the two DFT coefficients to the frequency was 
analytically derived and used in the algorithm. This is a new analytical result 
and forms the basis of the algorithm. 

PERFORMANCE ANALYSIS 

For the case of signal in additive noise, is a random variable. 
var[f.T.] = var[^]= vart^D,, )] 

DfoQ- 1 *'- 1 * 1 

In general, |a| and |0| are both Rician distributed random variables. However, 
under high signal to noise ratio, both |a| and |0| are essentially Gaussian 
distributed random variables. 

, which is a random variable, is found by the constraint, 

ff^will be perturbed by the noise component in D,,. Even though DJs 
constrained to be zero, the constraint and noise induce randomness in e,,. The 
noise perturbation in 
D. induces the perturbation in e„. 

The approach taken is the computation of the variance of Dfrom the point of 
view of the creation of D from noisy observations and then to find the 
corresponding perturbation of - e . 



For high signal to noise ratios, Dfo^Jmay be represented by the first three 
terms of the Taylor Series expansion about, p a and n p , which are the means 
of |a| and respectively. 
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Assuming e x =e 
A 

Ma =Mp =" 



srn( — ) 
V 2]ST 



5 and for high signal to noise ratio, 
I 

2 



a-l = al =-Ncr 



Then, 
E[D]=0, 

10 varTD]- 



It then follows that, 



Nsin 2 (— ) 

var[D] = 2^L, 

4(SNR) 



A 2 

15 where, SNR = Ar- 

<r 2 



£-£ = - tan*[D tan(— )] 
it 2N 



Then, for high signal to noise ratio, the normalized frequency error may be 
20 computed. The largest part of the probability density function of P is in the 
region of where the atan(x) « x . Therefore, 

al =vax(e-s) = 2N_ = 2N = 2N 2N 

4(SNR);r 2 



25 PERFORMANCE COMPARISON TO THE CRAMER-RAO LOWER BOUND 
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The Cramer-Rao Lower Bound on the variance of the frequency error of 
any unbiased frequency estimator is given by, 

a 2 = 6 

^ (2^) 2 N(N 2 -1)(SNR) 

The performance of the DFT based estimator may be compared to the Cramer- 
Rao Lower Bound. 

d Ni(N '- 1) ^ 2 O tm ' ( 4> 

GcRLB & 

For high SNR and iarge N, the performance frequency estimation variance is 
101og 10 (— ) = 0.063282577 dB .above the Cramer-Rao Lower Bound. 

_2 

Figure 5 shows in dB verses N, where N is the length of the FFT. 

JUSTIFICATION FOR ITERATION TOWARD A ZERO VALUE OF THE 

DISCRIMINANT 

The reason for the performance improvement of the proposed class of 
algorithms relative to prior algorithms is the first frequency interpolation allows 
the computation of two DFT coefficients, which are Y* DFT bin spacing above 
the first interpolated frequency and Vi DFT bin space below the first 
interpolated frequency. While the first interpolation may still have significant 
error, which is dependent on the relationship of the true frequency relative to 
the FFT coefficient frequencies, the error discriminant evaluated for the first 
interpolated frequency will have a value close to zero. The variance of the 
frequency error is relatively low in the region of small values of the frequency 
discriminant. Therefore, the second interpolated frequency will have small error 
variance. There is significant noise performance advantage in using the first 
interpolation to allow a low error variance second interpolation. The 
interpolation may be iterated, with diminishing Improvements of estimation 
accuracy, until convergence to a fixed point solution is obtained. Figure 6 
shows the variance of the normalised estimator frequency error estimate vs the 
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frequency error for the first interpolation. For the figure, N=64 and the signal to 
noise ratio is 6 dB. There is a very sharp reduction the rms error of the 
frequency estimator In the region of the frequency being close to the center 
frequency of the frequency discriminator. This indicates that tremendous 
5 improvement in performance obtained by Iteration. 

ALGORITHM SIMPLIFICIATONS RESULTING IN LARGE REDUCTIONS IN 

COMPUTATIONAL COMPLEXITY 



10 Simulation results have verified a single FFT and two iterations involving 

the computation of the discriminant function D(^^ 0 )andD(^,^ 1 )are sufficient to 

obtain RMS frequency error performance close to the computation of f . The 

first iteration moves the discriminant towards zero and decreases of | e-e x |. 

The estimate resulting from the second iteration therefore results in small error 
15 variance of s-e 2 . 

The algorithm has complexity of 0(N log 2 (N) +0(8N)=O[N log 2 (N)] 

The algorithm has the same order of complexity as the original FFT for 
performance, which is very close to the Cramer-Rao Lower Bound. 

The algorithm has the same order of complexity as the original FFT for 
20 performance, which is very close to the Cramer-Rao Lower Bound. 

Iteration with ^,(D) = D yield results very close to y(D) and saves 
considerable computational complexity. The fixed point solution for the two 
cases is identical. 

Since the algorithm comes very close to convergence with two iterations, 
25 two iterations are sufficient for most applications. The performance 
improvement with additional iterations Is small. 



OTHER FREQUENCY DISCRIMINANTS WITH THE SAME NOISE 

PERFORMANCE 

30 

There are a number of discriminants, which have the same performance, 
when used iteratively to obtain the fixed point solution, as the previously 
introduced discriminants. The noise performance is identical, for iteration, 
because the fixed point solution is identical. 
35 This class of discriminants includes functional forms, 
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and in particular, D=lL£LHgj! 

2|/?| 2 +|tf| 2 

and, 

D = Re ^+^^ ' wnere * denotes complex conjugate. 
And where Re[.] is the real part. 

SIMULATION 

Figures 7 to 12 show the rms frequency error performance of the 
algorithm vs SNR in dB, for N=2,4, 16,64,246, and 1024, respectively. Both the 
cases of one interpolation and two iterative interpolations are shown. The two 
interpolation case is essentially achieves the performance of the infinite 
interpolation case. 

In summary, a new, low computational complexity, class of algorithms, 
which interpolates the result of a FFT, has been presented for the precise 
estimation of frequency of a complex exponential function In additive Gaussian 
noise. The performance of the algorithm, above the threshold in additive 
Gaussian noise, is about 0.06 dB above the Cramer-Rao lower bound. The 
algorithm Is ideally suited to be utilized in a number of communications, signal 
processing and biomedical applications. The algorithm also has ideal 
characteristics for digital signal processor implementation. 

Industrial Application 

There are a large number of applications for this invention, including: 

- Rapid frequency initialisation of a phase lock loop for rapid signal 
acquisition; 

- Radar processing for precision radial velocity and radial acceleration target 
measurements; 

- Sonar processing for precision radial velocity and radial acceleration target 
measurements; 

- Satellite orbit determination; 

- Missile trajectory determination; 

- Ultra sound Imaging Doppler measurements for blood and other biological 
fluid velocity measurements; 
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- Ultra sound imaging for tomography processing involving Doppler shift 
measurements; 

- Coherent carrier tracking for coherent digital demodulators - large 
frequency acquisition range and rapid signal acquisition; 

5 - Noncoherent digital communication system frequency tracking - large 
frequency acquisition range and rapid signal acquisition; 

- Frequency estimation for electronic test equipment displays including 
frequency meters, oscilloscopes, spectrum analyzers and network 
analyzers; 

10 - Ultra low distortion, ultra high performance FM demodulator; and 

- Generalised software modules in Matlab and other commercial software 
packages. 

It will be appreciated by persons skilled in the art that numerous 
15 variations and/or modifications may be made to the invention as shown in the 
specific embodiments without departing from the spirit or scope of the invention 
as broadly described. The present embodiments are, therefore, to be 
considered in all respects as illustrative and not restrictive. 
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